USING THE SHERMAN-MORRISON- WOODBURY FORMULA TO SOLVE 
THE SYSTEM OF LINEAR EQUATIONS FROM THE STANDARD 
MULTIPLE SHOOTING METHOD FOR A LINEAR TWO POINT 
BOUNDARY- VALUE PROBLEM IS A BAD IDEA 

IVO HEDTKE 

Abstract. We use the standard multiple shooting method to solve a linear two point boundary- 
value problem. To ensure that the solution obtained by combining the partial solutions is 
continuous and satisfies the boundary conditions, we have to solve a system of linear equations. 
Our idea is to first solve a bidiagonal system related to the original system of linear equations, 
and then update it with the Sherman-Morrison- Woodbury formula. We study the feasibility, the 
numerical stability and the running time of this method. The results are: The method described 
above has the same stability problems like the well known Condensing method. The running 
time analysis shows that the new method is slower than the Condensing method. Therefore we 
recommend not to use the method described in this article. 



1. Introduction 

We solve the linear two point boundary-value problem 

££x(t) := x(t) - A(t)x(t) = r(t), t e [a, 6] 

3§x{t) := B a x(a) + B b x(b) = (3 

with the standard multiple shooting method, where x(t), r(t) : [a, b] — > R™, (3 £ R™, A(i) : [a, b] — > 
R" x ™ and B a , B b S R" xn . We divide the interval [a, b] with the shooting points 

a = r < 7i < . . . < TfYi — l < T m = b 

into m segments [V;, Tj+i]. We use the principle of superposition on each segment to find the 
solution 

Xj{t) = X{t;T 3 )c 3 +v{t;Tj), 
where Cj is a constant vector. X(t; Tj) is a fundamental system which fulfills the IVP 

jSfX(t; Tj) =0,te [Tj,T j+1 ], X^-Tj) = I. 

v{t\Tj) is an inhomogeneous solution of the ODE and fulfills 

3?v(t;Tj) = r(t), t e [Tj,Tj+l], v{Tj\Tj) = 0. 

The problem now consists in determining the vectors Cj in such a way, that 

(1) the function x(t) pieced together by the Xj(t) is continuous and 

(2) satisfies the boundary conditions. 

We define Xj := X(j-j-\-i',Tj) and Vj := v(rj+±;Tj). To satisfy the boundary conditions we focus 
on 38x{t) = j3: 

(1) B a co + B h X m -ic m _ 1 =P-B b v m -i. 

To ensure that x(t) is a continuous function we need 

x fc _i(r fe ) = XkCnO) k = l,...,m- 1, 
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which yields to the conditions 

(2) c fe - A"fe_ 1 c fe _ 1 = v k -i, k = l, ...,m-l. 

Now we collect equation ([TJ and the to — 1 equations @ in the following system of linear equations: 

(3) Mc = q, 
where we define Yj := —Xj and 
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It is known that M is regular if we assume that the 



Note that c, q € R mn and M e 
BVP has an unique solution. In this case 

(4) N := B a + B b X{b;a) 

is regular, too. (see Satz 8.1 (Theorem 8.1)]) 

2. The aim of this work 

There exists the well known method Condensing to solve the system (J3]) (see Section [6]) . Because 
of the special structure of M it is pretty obvious to try to find the solution in the following 
way: First solve the bidiagonal system from ((6]) and then update the solution with the Sherman- 
Morrison- Woodbury formula. In this paper we study the feasibility, the numerical stability and 
the running time of this method. 

3. The Sherman-Morrison- Woodbury formula 

Let A be a regular £ X I matrix and U and V be two I x p matrices. If I p + V T A^ 1 !! is regular, 
then 



(5) 

holds. 



(A + UV T y 1 = A' 1 - A~ X U{I P + V T A- 1 U)- 1 V T A' 1 . 



4. IS IT POSSIBLE TO USE THE SHERMAN-MORRISON-WOODBURY FORMULA TO SOLVE 

Mc = ql 

First, we have to split M into two matrices M = M, + U, where U can be written in the form 
U = UV T with U, V G W nnxn . For this we define 



U = 



0. 



,0,Bj 



and V T = [I n ,0,...,-L], 



where L := X n 1 • • • X,J_ . Therefore we have 



U = UV 1 = 







B„ 



and 



(6) 



where B := BhX 



b^-m-l 



M = M-U = 



B a L. 



B n L 



Y I 

Yx I 



Y m -2 I 

B 
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Now we have to check that M. is regular. Because of 

nm—2 
_ Q detYj, 

it follows that det A4 ^ iff dct B ^ 0, because the Yj are fundamental systems. But B = NL 
and TV and L are both regular. This follows from 



(7) 



X(b;a)=H m X 

A 7 — 1 



m-j- 



This shows that M. is regular. 

Finally we have to check that I n + V T M~ 1 U is regular. First we need an auxiliary result: 

Lemma. Given m regular n x n matrices Di . Then, the matrix 



A :-- 



D I n 

£>1 In 



D m -2 In 

Dm-1 



is regular and 



D7 1 



(D 2 D 1 D )- 1 
~{D 2 D 1 )- 1 



(_l)m- l(£) m _ 1 



Dr 



— {Drn-lDm-2) 1 



holds. 



Proof. It holds det A = J|"l 1 det£)j ^ 0. AA 1 = I mn and A X A = I mn can easily be 



verified. 



□ 



Now we go back to the matrix I n + V J M 1 U. With M.- 1 we denote the jth column of M. 1 
and we write M~[j for the n x n sub-matrix in the ith row and jth column of M^ 1 . With the 
lemma above and the new notation we get 

V T M- 1 U= [I n ,0,...,0 7 -L][M^\ ... I AC] [o,...,0,Bj 

= [M^ - LM£ I ... I - LM m U [0, • • • , 0, Bj 
= M^ l B a -LM m 1 m B a . 

With the special structure of M.^ 1 we can calculate the two sub-matrices and very 

easy: M^ m = B^ 1 and 



Mil = (-I)™- 1 isn^-i 

— 1 V 



(-ir- 2 Y r i ---Y-i 2 B- 1 



"771-2 



E _1 = LB' 



Now it follows that 

V T M- 1 U = M^Ba - LM m \ n B a = LB~ 1 B a - LB x B a = 0. 

The result above shows that I n + V T A4~ 1 U = I n is regular and we can use the Sherman- 
Morrison- Woodbury formula to solve ([3]) . 
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5. Solving Mc = q with the Sherman-Morrison- Woodbury formula 
With (J5J) the solution of © can now be expressed as 

c = M l q = {M+U)- X q = M^q - A^ _1 t7(J„ + V T M^U)' 1 V T M^q 
= MT x q - M- 1 UV T M- 1 q = MT x q - MT X UMT X q. 

This gives us an algorithm to solve (j3]): 

(1) Solve M£ = q. 

(2) Solve MC=U$. 

(3) Calculate c = £ — £. 

First we study the problem (1.) in detail. We have to solve 
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9m-l. 



Therefore we solve fi£ TO _j = <Z m _i and use recursion to find the other £ 

^i^j = 9j - Cj+i, 3 = m - 2, . . . ,0. 
We use the same method for our problem (2.). After we calculated 



u ... u 4 

m = 



B a B a L 

the resulting system of linear equations is 

V I 

Y x I 



£m-l. 







' Co ' 

Ci 

Cm-2 
Cm— 1. 



o 
o 



Again we first solve BC, m _ 1 = B a (£ — L£ m _ x ) and then solve the remaining systems of linear 
equations with recursion: 



Y jGj C, ,. j = m-2,. 
6. Condensing 



,0. 



We want to compare the new method above with the well known standard method from Stoer and 
Bulirsch. They solve ^ in the following way (see pQ or [3]): 

(1) Compute E := B a + B b X m - X ■ ■ ■ X and u := q rn _ 1 - B b X m -i(q m _ 2 + X m - 2 q m -3 + 
■ ■ ■ + X m -2 • • • Xiq Q ). 

(2) Solve Ec Q = u. 

(3) Compute the remaining Cj with recursion: c,-_|_i = q^ + XjCj. 

In the first step of our new algorithm from the section above we solve B£ m _ 1 = <? m _i- Notice 
that B = NL. But N = E holds. This follows directly from Q and 0. That means our 
new algorithm has the same stability problems like the Condensing method. See [jQ and [3] for a 
detailed discussion. 

Therefore we only analyse the number of flops used by the two algorithms to compare them. 
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Table 1. Running time analysis for the Condensing method. 




step 


description 


flops 




1 


Compute E and u. Because we compute the products of the Xj matrices in 
E we can use them to compute u, too. Therefore we need no extra product 








computations of matrices to compute u. 

• m — 1 matrix-matrix multiplications for E 

• one matrix addition for E 


(m - l)(2n 3 - 
n 2 


-n 2 ) 




• m — 1 matrix- vector products for u 

• m vector additions for u 


(m - l)(2n 2 - 
mn 


-n) 

J 


2 


Solve Eco = u. 


2/3n a 




3 


Compute the remaining Cj with recursion. 

• m — 1 matrix- vector products 

• m — 1 vector additions 


(m - l)(2n 2 - 
(m — l)n 


-n) 


E 


= 2mn i + 3mn 2 - 4/3n 3 - 2n 2 + n flops 






Table 2. Running time analysis of our new method. 


step 


description 


flops 




1 

1.1 
1.2 


Solve M£, = q. 
Solve B£ m _ a = q m _ x . 

• Compute T := L^ 1 = X m -2 ■ ■ ■ Xq. 

• Compute N ~ B a + B b X m -iT. 

• Solve Ns = q m _ 1 . 

. Compute £ m _i = iT 1 ^^ = L^N^q^ = Ts. 

Use recursion to find the other 


(m - 2)(2n 3 - 
4n 3 - n 2 
2/3n 3 
2n 2 - n 
(m-2)(2/3n 3 


n 2 ) 
+ n) 


2 

2.1 
2.2 


Solve Mi = W£. 

Solve BCm-i = B«(«o - ^m-l)- 
. Solve Tt = £ m _ r 

• Compute Ba(£ — *)• 

• Solve Ns = B„(£ - *)• 

• Compute C m _i = Ts. 

Use recursion to find the other £ ? . 


2/3n 3 

2n 2 

2/3n 3 

2n 2 - n 

(m-2)(2/3n 3 


) 


3 


Compute c = £ — £. 


mn 




E 


= lO/Smn 3 - mn 2 + mn - 2/3^ + 7n 2 - 4n flops 







7. Running time analysis 

We use LU- factorization to solve the systems of linear equations. We assume that this needs 2/3n 3 
flops for a n x n system. 

The running time of the Condensing method is analyzed in Table [TJ For a running time analysis 
of our new method see Table The result is: The Condensing method is faster than the new 
method described above. 

8. Conclusion 

We found a new algorithm to solve the system of linear equations from the boundary and 
continuity conditions with the Sherman-Morrison- Woodbury formula. This new method has the 
same stability problems like the Condensing method. Our new method is also slower than the 
Condensing method. Therefore it is not recommendable to use the Sherman-Morrison- Woodbury 
formula in this case. 
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